program ppCurrents

!ppCurrents.f90
!
! Last modified 13 August 2008
! Author: Cameron Bracken
!         Student Researcher
!         NOAA National Weather Service
!         300 Startare Drive
!         Eureka, CA 95501
!         cameron.bracken@noaa.gov
!
!Overview:
!   This program reads adcirc output and chops out a chunk based 
!   on a given bounding box.
!
!Longer Overview:
!   This program reads adcirc native velocity output (fort.64) 
!   and an adcirc grid file (fort.14) and ouptpust a subset of the 
!   velocity data and nodes at based on the bounding box in a file 
!   called 'domain'
!
!Inputs:
!   None, but the program looks for fort.14 and fort.64 in ../adcirc/
!   as well as looking for the domain file in chopped/
!
!Outputs:
!   NOTE: All outputs are written to chopped/
!   There are three output files:
!       fort.14.chop  - the node coordinates which fall in the the chopped region 
!       fort.64.chop  - identical to the fort.14 with the data falling outside 
!                       the region removed AND the time stamps removed 
!       fort.64.times - the timestamps which were removed from the fort.64
!
!


    implicit none
    double precision :: u,v,x,y,z
    integer, allocatable, dimension(:) :: insideNums
    double precision, dimension(2) :: xlim, ylim
    double precision :: ddum1,ddum2,ddum3
    integer :: idum1, idum2, nsets, nnodes
    integer :: i, j, k, n, eof, nInside, nTotalNodes, counter
    character(len=80) :: cdum

 
    ! Input files
    open(14,file="../adcirc/fort.14")
    open(64,file="../adcirc/fort.64")
    open(101,file="chopped/domain") 
        read(101,*)xlim(1),xlim(2),ylim(1),ylim(2)   
    close(101)

    ! Output files
    open(114,file="chopped/fort.14.chop")
    open(65,file="chopped/fort.64.chop")
    open(66,file="chopped/fort.64.times")

    
    !read coordinates from the fort.14
    
    read(14,*)cdum      !skip the first two lines
    read(14,*)nTotalNodes,idum1
    nInside = 0
    counter = 0 
    write(*,'(a)')"Searching for number of node in region... "
    do i = 1,nTotalNodes 
        if(mod(i,10000)==0)then
            write(*,'(1i3,a)')int(dble(i)/dble(nTotalNodes)*100d0),"% "
        end if 

        read(14,*,iostat=eof)idum1,x,y,z
        !write(*,*)idum1,x,y,z
        if(eof/=0)exit
        if(x>xlim(1).and.x<xlim(2).and.y>ylim(1).and.y<ylim(2))then
            nInside = nInside + 1
        end if 
    end do 
    write(*,'(a)')'Done.'
    write(*,*)'There are',nInside,'nodes in the region.'
    
        !now we know how many nodes in ther region there are, we can allocate
    allocate(insideNums(nInside))
        
        !go back through and get the actual coordinates
    rewind(14)
    read(14,*)cdum      !skip the first two lines
    read(14,*)idum1,nTotalNodes

    n = 1
    insideNums = 0
    write(*,'(a)')"Getting node numbers in region... "
    do i = 1, nTotalNodes
        
        if(mod(i,10000)==0)then
            write(*,'(1i3,a)')int(dble(i)/dble(nTotalNodes)*100d0),"% "
        end if 

        read(14,*,iostat=eof)idum1,x,y,z

        if(x>xlim(1).and.x<xlim(2).and.y>ylim(1).and.y<ylim(2))then
            write(114,*)x,y,z
            insideNums(n) = i
            n = n + 1
        end if 

    end do
    write(*,'(a)')'Done.'
    !write(*,*)'Wrote nodes to nodes.'
    
    !Now read the fort.64 file
    
    write(*,"(a)")'Writing chopped region... '

    read(64,*)cdum      !skip the first line
    write(65,*)'Skip this line'

    read(64,*)nsets,ddum2,ddum1,idum1,idum2
    write(65,*)nsets,nInside,ddum1,idum1,idum2

!write(*,*)insideNums

    do i = 1, nsets
        
        if(mod(i,10)==0)then
            write(*,'(1i3,a)')int(dble(i)/dble(nsets)*100d0),"% "
        end if 

        n = 1
        read(64,*)ddum1,idum1
        !write(*,*)ddum1,idum1
        write(66,*)ddum1,idum1

        do j = 1, nTotalNodes
            read(64,*)idum1,u,v

            if(insideNums(n)==j)then
                write(65,"(1i10,2g20.10)")n,u,v
                n = n + 1
            end if 

        end do 
    end do
    write(*,'(a)')'Done.'
    write(*,*)'Wrote chopped region to fort.64.chop and times to fort.64.times.'

    close(14)
    close(114)
    close(64) 
    close(65)
    close(66)

    stop
end program ppCurrents
